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Abstract 

When using images to locate objects, tiiere is the problem of correcting for dis- 
tortion and misalignment in the images. An elegant way of solving this problem 
is to generate an error correcting function that maps points in an image to their 
corrected locations. We generate such a function by fitting a polynomial to a set of 
sample points. The objective is to identify a polynomial that passes "sufficiently 
close" to these points with "good" approximation of intermediate points. In the 
past, it has been difficult to achieve good global polynomial approximation using 
only sample points. We report on the development of a global polynomial approx- 
imation algorithm for solving this problem. 

Key Words: Polynomial approximation, interpolation, image rectification. 



1 Introduction 

The problem that is addressed here occurred in the context of the development of a 
simple, low-cost robotic exhibit for demonstrating the concept of intelligent robotics to 
the general public. Intelligent robotics deals with the use of sensors to enhance a robots 
performance in an uncertain environment. For the exhibit, we implemented a visually- 
guided pick-and-place robot. The robot uses an image to determine the location of 
objects placed arbitrarily on a flat surface and demonstrates success in locating the 
objects by manipulating them. The exhibit consists of a robotic arm that is fixed in 
front of a flat work surface. Two pedestals are placed anywhere within reach of the 
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arm. One pedestal is blue and the other is green. A blue ball is placed on the blue 
pedestal and the robot must pick up the ball and place it on the green pedestal. An 
ordinary webcam is placed in a frame above the work surface and is used to determine 
the location of the pedestals so that the arm can be guided accordingly. Fig. [Tf a) shows 
the camera's view of the workspace. Once the exhibit has been set up, camera, robot 
and workspace are all fixed relative to one another 




(a) unrectified image (b) rectified image 



Figure 1: Rectification of an image of a test pattern placed over the robot's workspace taken 
from a severely misaligned camera. A grid has been superimposed on the images to aid in 
comparing the horizontal and vertical afignment of image features. Although the original images 
are colour coded, the images shown are monochrome. If the upper left dot and the upper right 
dot are ignored, the rest of the dots in (b) sample the area of the workspace that is within reach 
of the robot. 

There is therefore the problem of determining the location of the colour-coded 
pedestals based on their image. This problem is compounded by the fact that the we- 
bcam produces significant image distortion, and by the fact that the position and ori- 
entation of the camera relative to the work surface may not be exactly the same every 
time the exhibit is set up. A similar problem occurs when any automated mechanism 
is taken apart for maintenance. The mechanism usually must be recalibrated when it is 
reassembled. In this sense, we are addressing the problem of easy recalibration of the 
vision component of our robotic exhibit. 

In processing the image, there is a need to rectify significant image distortion 
caused by the optical properties of the camera and by errors in positioning the cam- 
era during set up of the exhibit. We solve this problem by determining a mapping from 
the pixel locations of image points to the physical locations of corresponding source 
points. The mapping has to be determined empirically in order to recalibrate the vision 
system each time the apparatus is set up. The pixel positions of the images of key 
points in a test pattern are mapped to the known locations of these key points. This 
partial mapping is then used to approximate a mapping for the entire image. 

In |T|, Brown surveyed and classified several established methods for determining 
the mapping from pixel position to source location for a camera. Using Brown's clas- 
sification, the errors caused by distortion and misaUgnment are static in the sense that 
they do not change from image to image taken with the same camera in the same posi- 
tion. Static distortions can be rectified in a one-time-only setup process via calibration 
techniques. Our method may be regarded as a calibration technique. 

Brown classifies the distortion due to camera properties as internal and the mis- 
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alignment as external. Our approach does not require the use of any sort of model of 
the characteristics of the camera or the geometry of how an image is captured. We treat 
the mechanism by which source points are captured as image points as a "black box", 
the inner workings of which is not modelled. We therefore handle these internal and 
external distortions without having to distinguish between them. 

Furthermore, we are concerned with geometric distortions as distinct from photo- 
metric distortions. The effect of photometric distortions are minimised by using pri- 
mary colours (red, green, blue) for the key features in the scene and either black or 
white for the other features. Features can therefore be identified by colour without re- 
gard for intensity values. Colour distortions are handled by filtering out colours below 
a fixed saturation value using a median filter to remove speckles and discretising the 
remaining colours to fully saturated red, green or blue Q. Visual artifacts are then 
located by looking for regions within the image with the appropriate colour. 

The problem is thus reduced to one of finding an adequate approximation of a total 
mapping from pixel positions to object locations using a set of sample points. Three 
popular approaches to solving this problem are: the use of an artificial neural network; 
the use of a piecewise interpolation technique to fit curves of a chosen form in between 
the sample points; the use of polynomial interpolation to fit a polynomial to the sample 
points. 

Artificial neural networks are popular among engineers for their ability to gener- 
aUse from sample data. We find it difficult to argue on practical grounds against this 
approach. However, for us, the main difficulty with this method is its inability to yield 
a representation of the approximating function that is independent from that of the neu- 
ral network itself. However, recent advances in the use of algebraic training of neural 
networks to produce closed form analytic solutions are addressing this problem ||3]|4|. 

Piecewise interpolation is popular among the data visualisation community for its 
ability to handle localised distortions 11]. Given sample points on a planar surface, 
the domain must be divided into polygonal segments with sample domain values as 
vertices. A smooth surface segment is then fitted over each polygon such that the 
sample points associated with the vertices lie on the surface, and the edges associated 
with adjacent segments meet in a prescribed way. There is usually more than one way 
to segment the domain space, and different segmentations may yield very different 
interpolations. 

Global (as opposed to piecewise) polynomial interpolation is a theoretical possibil- 
ity that has proven elusive in practice. The Stone- Weierstrass Theorem 15] [6l estab- 
lishes the existence of a polynomial approximation for every real-valued continuous 
function defined on a closed interval. Therefore, as long as the sample data does not 
admit the existence of discontinuities, we should be able to fit a polynomial to the 
points with an arbitrary degree of precision. The problem (exemplified by Runge's 
function fS", T]) is that successive interpolations do not necessarily converge as more 
data points are added to the set of sample points (c.f. Fabers Theorem ||5i|7J). We will 
refer to this problem later as the convergence problem. 

Our solution is similar to global polynomial interpolation. However, the fact that 
we are dealing with noisy data allows us to relax the requirement for an exact fit and to 
focus instead on the suggested shape of the function. It is because of this relaxed focus 
on an exact fit that we refer to our problem as an approximation problem rather than an 
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interpolation problem. By focusing on the suggested shape of the function, we directly 
address the convergence problem. 

In this paper, we present the generic function approximation algorithm and use 
image rectification as an example of its use. In section|2l we develop the theory behind 
the approximation algorithm in the univariate case. Section [3] discusses the bivariate 
version of the approximation algorithm that is used to solve our robot vision problem. 
Section|4]summarizes our findings. 

2 Polynomial approximation of functions 

In the univariate case, the polynomial approximation problem can be defined as fol- 
lows: 

Definition 1 (The univariate polynomial approximation problem) Find a polynomial 
P{x) that fits a set of m sample points {xi, Tji); where i ~ 1, 2, . . . , to and the Xi are 
distinct; such that: 

max \yi — 'P{xi)\ < e, for e > 

i—l 

That is, we are trying to identify a polynomial, P(a;), that fits a finite set of points to a 
degree of precision determined by e. The condition that the Xi must be distinct ensures 
that there are no discontinuities in the target function. Otherwise, a solution is not 
assured. 

The Weierstrass approximation theorem |I3] |5] |3 assures us that this is a solvable 
problem for e > 0: 

Theorem 1 (Weierstrass) IfF is any continuous function on the finite closed interval 
[a,b], then for every e > there exists a polynomial P„(a;) of degree n (where n 
depends on e) such that: 

max |F(a;) - P„(a;)| < e 

Our finite set of points is not generated from a known function that we are trying 
to approximate, so we will be looking instead for an approximation that generates 
intermediate points that are "good" in some generic sense. We require the algorithm 
to yield an approximating polynomial that does not "curve" unnecessarily between 
sample points. In our algorithm, we try to conservatively fit a shape to the sample points 
using a linear combination of Chebyshev polynomials. A preference for fitting lower 
order Chebyshev polynomials reflects a preference for simple shapes. This emphasis 
on simple shapes addresses the convergence problem by only adding detail when it 
results in a better fit. 

Section 12.11 presents the basis for our algorithm in its purest form and explains 
why it is inadequate. We refer to this algorithm as the Cartesian vector based (CVB) 
interpolation. This algorithm is of theoretical interest only, since it produces results 
similar to Lagrange interpolation, with the same convergence problems illustrated by 
Runge and covered by Fabers theorem [5,7]. It sets the scene for presentation of the 
modified version of the algorithm. In section 12.41 we present the modified algorithm 
that gives better results on the same data in terms of capturing the suggested shape of 
the curve. We refer to this second algorithm as the CVB approximation algorithm. 
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2.1 Casting the problem in Cartesian vector space 

We address the approximation problem by finding a linear combination of the first n 
Chebyshev polynomials of the first kind that exhibits the desired properties of P(a:;). 
Thus: 

n-l 

P(xO - ^a,T,(x) (1) 

i=0 

where Ti{x) is defined for a; e [—1, 1] as 

• To (a;) = 1, 

• Ti(a;) = X, 

• Ti{x) = 2xTi_i{x) - Ti_2{x), fori > 1. 

The problem is to find a set of coefficients, at, i = 0, . . . ,n — 1 that establish a fit. 
Although other sets of basis polynomials exist, Chebyshev polynomials reputedly yield 
good interpolation results. In what follows, we will assume without loss of generality 
that X falls within the interval [—1,1]. 

The fact that we are using a finite set of sample points allows us to reformulate the 
problem as one involving vectors in a Cartesian m-space; where m is the number of 
sample points. For this purpose the problem is restated as follows: 

Deiuiition 2 (The univariate Cartesian vector based (CVB) polynomial approximation problem) 

Given 

• a set of points, {xi,yi) G M?, i = 1, . . . ,m; 

• a set of m-dimensional Cartesian vectors, tj, j = 0, . . . ,n — 1, such that the 
i'th component ofrj is equal to Tj{xi),for i = 1, . . . ,m; 

• a Cartesian vector, 7 such that the i'th component 0/7 is equal to yi, for i = 
1, . . . ,to; 

• a real value e fe > Oj; 

find values for a set of scalar quantities, aj, such that: 

En — 1 

• (5 = 7 — p; 

• l^'l < ^for each component, 6^, of 6 (i = 1,. . . , m). 
For example, the representation for the expression 

x^ + 2x + l 

as a Unear combination of Chebyshev polynomials is: 

1.5To(;c) + 2Ti(x)+ 0.5X2(0;) 
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Table 1 : The CVB representation of a curve fit. 
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Table[T]shows the situation for this function for three evenly spaced sample points. 
In this example, each Chebyshev polynomial yields a 3-dimensional Cartesian vector. 
A linear combination of these vectors yields a solution. We will refer to the vectors 
generated by the polynomial terms as term vectors. 

If the term vectors tj are orthogonal, they may be taken as basis vectors and solu- 
tion values for the coefficients, aj, may be directly obtained from the projection of 7 
on each tj, respectively. 

However, orthogonality does not hold in general, as illustrated by the set of points 
in Table [T] An obvious solution to this dilemma is to first identify an orthogonal set of 
vectors corresponding to the term vectors and use this set to determine a solution. We 
developed an algorithm to do just this. We refer to it as the CVB interpolation. 

2.2 The CVB interpolation 

This algorithm computes an orthogonal set of vectors, {oq, oi, . . . , Om_i}, from the 
first m linearly independent term vectors, derives the projection of 7 on each of these 
vectors and modifies the coefficients of the term vectors accordingly. 

The orthogonalisation of the term vectors is based on the following theorem on 
which the well-known Gram-Schmidt orthogonalisation process is based: 

Theorem 2 Let {vi, V2, . . . , v„} foe an orthononnal set of vectors in a vector space, 
V. Then for any vector w G the vector: 

n 

O = W - • Vi)Vi 

i=l 

is orthogonal to each v^, i = 1, . . . , n. 

This theorem can easily be generalised to apply to orthogonal (as distinct from or- 
thonormal) sets as follows: 

Theorem 3 Let { vi , V2 , . . . , v„ } foe an orthogonal set of vectors in a vector space, V. 
Then for any vector w ^ V, the vector: 

E(W • Vi) 
7 

2—1 ^ ^ 

is orthogonal to each Vi, i — I, . . . ,n. 
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The additional factor in the summation normalises each v^. 

In our algorithm, we identify an orthogonal set of vectors, {oq, Oi, . . . , o„_i} and 
a set of scalars pj^k, such that: 

oo = To (2) 

_ (Tj ■ Ok) 

Pj,k - 7 r (3) 

°J = 'Tj ~^Pj,kOk,foij ^l,...,n-l (4) 
fe=o 

We refer to each Oj as an orthogonal component of the corresponding t,, for j = 
0, ...,n- 1. 

For the CVB interpolation, we need to express Oj as a function of Tfc. Using (|4|i we 
write: 

3 

k=0 

k=Q 

k=Q 1=0 
1=0 k=l 



A:=0 /=/c 



From (|5]l we get: 



Ijj = 1 (6) 



Fig. |2] shows the CVB interpolation in pseudocode. Table |2] shows a trace of the 
CVB interpolation on the problem depicted in Table [T] 

2.3 Difficulties with polynomial interpolation 

The following types of data serve to illustrate the difficulties encountered with global 
polynomial interpolation lH: 

• "Humped and flat data" that suggest a flat curve in some regions and not in others 
(Fig.EJ. 
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Computs th© o]rtho<^on3l cornponsnts, Oj and ths scslairs Qj k 


3 . 


For j = 0,...,n—l do 




(a) Compute the error vector, 5. 




(b) Set Oinc to 




( c ) For k — Q, . . . , j do 




Set Ofc to Ofe + Oinc ■ Qj^fe . 




Endf or 




Endf or 


4 . 


Return the set of coefficients, aj . 



Figure 2; The CVB interpolation (see section |22] |. Computational details have been omitted 
that deal with avoidance of representational and computational error. 



Table 2: A trace of the CVB interpolation on the problem depicted in Table[T] 





ai 


a2 


l|5|| 


P(^i) 






1.666666667 


0.0 


0.0 


2.943920289 


1.666666667 


1.666666667 


1.666666667 


1.666666667 


2.0 


0.0 


0.816496581 


-0.333333333 


1.666666667 


3.666666667 


1.5 


2.0 


0.5 


0.0 


0.0 


1.0 


4.0 



• "Noisy straight line data" with y values given at unevenly spaced x values (Fig. 

in. 

• Data from functions that exhibit the convergence problem, the classical example 
of which is the Runge function (Fig. |5]). Faber's theorem establishes that the 
Runge function is not the only function that exhibits this problem. 

These types of data all exhibit a common problem: there is significant deviation of 
the fitted polynomial from the shape of the curve suggested by the data points. With 
regard to the convergence problem, progressively adding more data points does not 
lead to progressively better agreement between the fitted polynomial and the suggested 
shape of the curve (c.f. Faber's Theorem). The problem of large deviations from 
ground truth at the intermediate points is usually addressed in one of two ways: either 
additional information is supplied as data, or additional constraints are imposed as part 
of the method of solution. 

Hermite interpolation is an attempt to improve the quality of the interpolation by 
using additional information (first derivatives at the data points), but suffers from the 
fact that the additional information is local to the data points and so does not sufficiently 
constrain what happens at the intermediate points. What seems to be needed are global 
constraints that affect the quality of the solution over the entire interval. 
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(a) interpolation on 9 points (b) approximation on data from (a) 

Figure 3: Interpolation and approximation of data that is flat in some places and humped in 
others. 




(a) interpolation on 9 points (b) approximation on data from (a) 



Figure 4: Interpolation and approximation of noisy straight line data with y values given at 
unevenly spaced x values. 

The difficulty with applying global constraints in interpolation is the requirement 
that a perfect fit be achieved at the sample points. This requirement gives the sample 
points a special status relative to the intermediate points. The fact that a constraint must 
respect this special status would render it less than global, as the constraint would have 
to behave differently in the vicinity of the sample points. If a perfect fit is not required, 
then global approximation may be used. Theory states that polynomial approximations 
exist that follow closely any curve that can be described by a continuous real-valued 
function on a closed interval. Global constraints can be more easily applied to selecting 
the best approximation than to finding a perfect fit. 
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(c) approximation on data from (a) (d) approximation on data from (b) 



Figure 5: Performance on thie convergence problem. Data sets were taken from the Runge 
function. Comparison of (a) with (b) illustrates that interpolating with more data points does not 
always give better results. Comparison of (c) and (d) illustrate convergence on the same data 
using the CVB approximation algorithm presented in section l2!4l 

2.4 The CVB approximating algorithm 

In this section, we present our CVB approximating algorithm. This algorithm gives 
better convergence than the CVB interpolation. This improvement was achieved by 
exploiting the observation that the shape of a function that fits the error vector, S (see 
Definition |2|l, is captured in part by the direction of S, together with the fact that the 
direction of a Cartesian vector is invariant under scalar multiplication. We define the 
shape of a function as follows: 

Deflnition 3 (The shape of a function) Consider two continuous real-valued functions, 
Fi and F2, defined on the same closed interval, I. Fi and F2 are of the same shape if 
there exists a finite non-zero real scale factor, s, such that Fi ~ SF2 

Thus, a particular shape is denoted by the set of all the functions that are of that shape. 
We will only explore in this paper the properties of this concept of shape that are 
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relevant to our exposition. In the first instance, we will confine our discussion to real- 
valued functions defined on the interval / = [—1,1]. 

Definition 4 (Value ratios) If two functions, Fi and ¥2, are of the same shape, then 
for any two distinct values, Xi,X2 G such that Fi(x2) 7^ f and, by extention, 
F2 {^2 ) ^ 0), it is the case that pj|||^^j = That is, the ratio of Fi {xi ) to Fi (x2 ) 

is the same as the ratio of F2{xi) to F2{x2)- 

It is this invariance of value ratios across functions of the same shape that defines our 
concept of shape. 

Now consider a Chebyshev polynomial, T, and its associated term vector, r, de- 
fined on a set of sample points. The ratios of the components of r to one another are 
value ratios of T. But these ratios also define the direction of the Cartesian vector, r. 
So that the shape of T determines the direction of r. However, since the finite dimen- 
sional term vector r cannot capture the infinity of value ratios that specify the shape of 
T, the direction of r only partly captures the shape of T. 

If we wish to compare the shape of a Chebyshev polynomial, T, with the shape 
suggested by 6, we can compare the direction of the corresponding term vector, r, with 
the direction of 6. In particular, we look at the projection of (5 on r and express the 
projection as a scalar multiple of r. The scalar factor thus identified is used to adjust 
the coefficient of T in equation ([T]i. 

The main difference between the CVB approximating algorithm and the CVB in- 
terpolation is an emphasis on fitting the shape of the term vectors as distinct from an 
orthogonal component thereof. However, since term vectors are not orthogonal in gen- 
eral, more than one approximation may be possible. As a means of selecting a preferred 
solution, we adopt an heuristic based on fitting lower order terms first. 

Like the CVB interpolation, the CVB approximation algorithm works by finding 
projections of the error vector on a vector space identified by a finite number of term 
vectors. It gives preference to using the term vectors associated with lower order 
Chebyshev terms. The process is as follows (see Fig. |6]for the algorithm in pseu- 
docode): 

• The first projection to be eliminated is the projection of 5 on tq. Whenever the 
projection of J on a term vector, r, has been eliminated, r will be referred to as 
having been visited. 

• If Tj has been visited, then Tj_i must be re-visited. This rule is applied recur- 
sively to revisits until tq is revisited. 

• After Tj has been visited and all consequent revisits have taken place, then tj^i 
is visited. 

• The process terminates when the required degree of accuracy of fit has been 
obtained. 

Revisits are necessary in order to maintain the fit of lower order terms. More preferred 
terms are revisited after lesser preferred ones in order to minimise the effect of revisits 
on their fit. Thus the most preferred term is revisited last. 
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2. Compute the error vector, S. 

3. Set j to 0. 

4 . While j < n and max(|5'|) > e for z = 1, • • • , m do 

T • -6 

(a) Set a,- to a, -| — - — . 

Tj-Tj 

(b) Compute the error vector, 5. 

(c) For fc = (j- l),0--2),...,0 do 

Compute the error vector, 5. 
Set Ofc to Ofe + . 
Compute the error vector, 5. 
Endf or 

(d) Set j to j + 1 
EndWhile 

5. Return the set of coefficients, aj . 



Figure 6: The CVB approximation algoritiim (see section l2!4t . Computational details have 
been omitted that deal with avoidance of representational and computational error. 



2.4.1 Convergence of the CVB approximation algorithm 

To prove convergence, we consider the general situation in which we have visited the 
first p term vectors and still have a non-zero residual error vector We note that we can 
eliminate this error if we can approximate the error vector using a linear combination 
of the Chebyshev polynomials that have not yet been visited. We use the real-valued 
version of the Stone- Weierstrass theorem f6l to prove that this is possible. 

Theorem 4 () Let X be a compact set and let C{X) denote the space of continuous 
real-valued functions defined on X. Assume that A is a subalgebra of C{X). Then A 
is dense in C{X) in the uniform norm iff A separates points and for each x ^ X there 
exists anf £ A satisfying f(x) ^ 0. 

Proof 

To prove that the polynomial vector space P with basis set {T„ \ n > p} is dense in 
C{X), where X is the interval [-1,1], prove that: 

1 . X is compact. 

2. P is a subalgebra of C(X). 

3. P separates points (for any distinct points, Xi,X2 G X , there is a T„, n > p s.t. 

T„(a;i) Tn{x2). 

4. For each x £ X there is a T„, n > p s.t. T„(a;) ^ 0. 
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Taking each condition in turn: 

1. X is a closed and bounded interval on the real number line and is therefore 
compact. 

2. P is a subalgebra of C{X) since P is closed for the usual multiplication and ad- 
dition of functions, and for scalar multiplication by real values. This is sufficient 
to ensure that the defining properties of an algebra hold for the subspace, P as 
they do for C{X). 

3. Let xi,X2 G [—1, 1] such that xi ^ X2 and T„(xi) = T„(a;2) . Given that 
Tn(x) = cos(ncos^^(a;)), let = cos(6'i) and X2 = cos(6'2) with 6*1 > 62 and 
6*1, 6*2 € [0, tt]. If T„(a;i) = T„(a:;2) , we get cos(n^?i) = cos(n02)- This implies 
that either n^i (mod 27r) = n92 (mod 27r) or —{ndi (mod 27r) = 27r — (71^2 
(mod 27r)), which gives 

2ttN 

Oi — 02 = for some integer, N (8) 

n 

Since 61,62 G [0, tt] and di > 62, it follows from © that: 

0<^<i (9, 

n 2 

Now, T„+i(xi) = cos((n + l)6i). So that from 

{n + 1)61 = 271 (^N + + {n + 1)62 (10) 

From (|9]l and ( fTOl i we can conclude that cos((n + 7^ cos((n + 1)62). So 

that Tn+i separates the points. 

4. The zeroes of T„ are given by (j"' = cos (^-^^^jn^) ^ 1, . . . , 71. So that 

T„ (cj"^) - 0. Let ^5") - i22_i)ZL so that c]"' - cos (^f If cj"^ is also 

a zero of T„+i then T„+i (^cj"') = cos (^(n + l)6lj"'^ = 0. This implies that 

^j"-* is a multiple of tt. However, by definition, this is not the case. So cj"-* 
cannot also be a zero of T„-|_i. Therefore, for any x G [—1, 1], if Tn{x) = 0, 
then Tn+i(x) and the result follows. 



Since we are solving the problem in the Cartesian m-space identified by the m 
sample points, we must prove that the subspace identified by the set of r„, where n > p, 
contains term vectors a linear combination of which will yield a vector arbitrarily close 
to the error vector. 
Proof 

Since the term vectors are derived from the Chebyshev polynomials and the subspace 
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P is dense in C{X), it follows that there is a linear combination of Chebyshev polyno- 
mials in this subspace that is arbitrarily close to a function passing through the points 
of the error vector. This will identify a Unear combination of term vectors that is arbi- 
trarily close to the error vector There will therefore always be at least one more term 
vector that reduces the magnitude of the error vector. 



As for the rate of convergence, it is a well-known result in approximation theory 
that if TO is the number of sample points, good results can be obtained by placing your 
sample points at the zeroes of Tm+i- This yields to orthogonal term vectors. For 
these points, the CVB approximation algorithm yields the same solution as the CVB 
interpolation algorithm. Furthermore, the quality of the approximation is consistent 
with the shape fitting properties of the CVB approximation algorithm. 

Since we were using noisy data, we did not exploit this result in our application. 
Instead, we used heavy sampling of the particular area of interest (the region of the 
workspace that the robot can actually reach) with minimal sampling of points at the 
Umits of the x and y intervals (i.e. the corners of the rectangular workspace). In general, 
because the norm of the residual error vector decreases as new terms are fitted, and in 
most cases only a component of the error vector is removed with each term visited, 
convergence slows down as the algorithm progresses. 

3 The bivariate CVB approximation algorithm 

We use a bivariate version of the CVB approximation algorithm to map image pixels 
onto object locations. The bivariate approximation problem is defined as follows: 

Definition 5 (Tlie bivariate polynomial approximation problem) Find a polynomial 
P{x, y) that fits a set ofm sample points {xi, yi, Zi); where i = 1, 2, . . . , m and the 
ordered pairs {xi,yi) are distinct; such that: 

max \zi — 'P{xi,yi)\ < e, for e > 

Thus: 

n — 1 n — 1 

P(2;,2/) = EE«'.^Ti(^)T,(2/) (11) 
i=o i=o 

where Tfc is the (k-Hl)th Chebyshev polynomial of the first kind; x,y & [—1, 1]; ^ = 
if i + j > n (i.e. the coefficients form a triangular array). As defined, P(a;, y) is a 
bivariate polynomial of degree n—1. 

Casting this problem in Cartesian vector space, we define the bivariate CVB poly- 
nomial approximation problem as follows: 

Definition 6 (The bivariate CVB polynomial approximation problem) Given 
• a set of points, {xi, yi, Zi) G K^, i = 1, . . . ,m; 
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• a set of m-dimensional Cartesian vectors, Tj^k, for j + k < n and j, k = 
0, ... ,71 — 1, such that the i'th component of Tj^k is equal to j{xi)Tk{yi), 
for i = 1, . . . , m; 

• a Cartesian vector, 7 such that the i 'th component of 7 is equal to Zi, for i = 
l,...,m; 

• a real value e fe > Oj; 

find values for a set of scalar quantities, aj^k, such that 

• P '^],k'''j,k for j + k < n and j, fc = 0, . . . , n — 1; 

• S — J — p; 

• I'^N — s for each component, S^, of 5 (i ~ 1, . . . , m). 

So that the bivariate CVB polynomial approximation problem involves a search for 
a linear combination of m-dimensional Cartesian vectors; where m is the number of 
sample points. Once the problem has been cast as a CVB problem, the coefficients of 
the solution are found through a progressive reduction of the residual error vector as 
for the univariate CVB polynomial approximation problem. 

In the univariate version of the CVB approximation algorithm, terms were visited in 
order of increasing degree of the corresponding Chebyshev polynomial. In the bivariate 
case, it is not immediately obvious how to define the order of preference of terms, since 
several terms can be of the same degree. The following preferencing rules produce 
acceptable results for our application: 

Let {i,j) denote the term with coefficient ai^: 

• The term (a, b) is visited before term {c, d) if a + b < c + d. That is, i + j is 
used as a primary measure of the "preference" for fitting one term over another. 

• If a + b = c + d, then (a, 6) is visited before term (c, d) if min(a, b) < min(c, d). 
That is, min(i, j) is used as a secondary heuristic preference metric. 

• Ifa + b^c+d and min(a, b) = min(c, d), then (a. b) is visited before term 
(c, d) if a < c. This is an arbitrary ordering rule to sequence terms of equal 
preference. 

• After term is visited, only terms (a, b) where a < i and b < j are revisited 
and revisits take place in reverse order to visits. 

Informally, when visiting, these rules give preference to a term that is of lower degree 
in the first instance, or if of equal degree, has a component that is of lower degree 
than any component of a lesser preferred term. When revisiting, the same preferences 
are observed, but revisits are restricted to terms with components no greater than the 
corresponding components of the last visited term. So that in general, preference is 
given to lower order terms. Proof of convergence is similar to that for the univariate 
case given in section l2".4.1l 
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3.1 Using the bivariate CVB approximation algorithm 

Fig. [T](a) shows a deliberately severe misalignment of the camera. The image of the 
workspace appears to be rotated in a clockwise direction. The camera also produces a 
pincushion distortion that is reportedly imperceptible to most observers but is signifi- 
cant with respect to locating objects on the workspace. If uncorrected, distortion results 
in an error of as much as eight millimetres in the location of an object (corresponding 
to a misplacement of four pixels in the location of an image point). 

In comparison, Fig. [r|(b) shows the results of using the bivariate CVB approxi- 
mation algorithm to rectify the image shown in Fig. [T](a). The rectified image shows 
a rectangular border to the workspace with edges that are better aligned horizontally 
and vertically, as can be seen by comparing their alignment with the superimposed 
grid. In practice, the misalignment will also be imperceptible, but both distortion and 
misalignment are sufficient to cause errors in locating objects on the workspace. 

No attempt was made to ensure that the key points in the test pattern were evenly 
spaced, or placed according to the zeros of a Chebyshev polynomial. Instead, all but 
two of the key point^ were used to sample the region of the workspace that is within 
reach of the robot. The two extra points were used primarily to fix the corners of the 
workspace for illustrative purposes. 

An even spacing of the key points in the test pattern is of questionable utility since 
this does not guarantee an even spacing of their images due to distortion and misalign- 
ment. Fig. m illustrates the stability of approximations as compared to interpolations 
when unevenly spaced points are used. This characteristic of approximations was ex- 
ploited here. 

As for using the zeroes of a Chebyshev polynomial, there is the question of which 
Chebyshev polynomial to use, since, for the CVB approximation algorithm, it is not 
predetermined how many Chebyshev terms will be included in the approximation. For 
the rectification depicted in Fig. [T] the CVB approximation algorithm yielded a so- 
lution with 27 bivariate Chebyshev terms. This yields a polynomial of degree 7. A 
bivariate polynomial of degree n can have up to ^{n + l)(n + 2) terms. Given that 
twenty sample points were used and a bivariate polynomial of degree 20^ — 400 can 
have up to 80601 terms, this represents a significant saving in computational overhead. 

4 Summary and conclusions 

We have developed what we believe to be a novel algorithm for global polynomial 
approximation. We call this algorithm the Cartesian Vector Based approximation algo- 
rithm, or CVB approximation algorithm. This algorithm has several desirable features: 

• It does not require one to fix the degree of the approximating polynomial ahead 
of time. The algorithm is capable of progressively adding polynomial terms until 
the required precision of fit is achieved (or some specified limit on resources is 
reached). 

' The key points are the centres of the solid circles in the test pattern shown in Fig. [T] 
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• Since the algorithm progressively converges on a "closest fitting function", an 
approximate solution is available after each iteration. The longer the algorithm 
runs, the better the fit. 

• The algorithm yields better results on the type of data that usually presents diffi- 
culties for global polynomial interpolation. 

It has been proven that the algorithm will yield an approximation that is within e of a 
perfect fit, where e > and the uniform norm is used as the distance metric. The rate 
of convergence depends on the choice of data points. 

Our focus has been on presenting the CVB approximation algorithm and we use 
image rectification to illustrate its use. We show how global polynomial approximation 
can be used to calibrate the vision component of a visually guided pick-and-place robot. 
Calibration problems are sufficiently prevalent within the field of robotics to render this 
example relevant. 

There are, however, other types of mappings within the field of robotics to which 
the algorithm may not be directly applicable. For instance, in time series analysis a 
study is made of a time-varying information-carrying signal for the purpose of predict- 
ing its future behaviour. Since time series analysis attempts to predict future events, 
there is an emphasis on extrapolation. Our algorithm is based on interpolation rather 
than extrapolation. It is an open question whether it can be adapted for time series 
analysis. 

Systems identification is another area that uses function approximation [9|. For 
instance, ARMAX/NARMAX models are parametrised models of systems consisting 
of time-varying input values and output values, where the assumption is that the output 
values depend in some way on current and past input values. System identification 
involves finding an instantiation of the parameters that yields a predictor of system be- 
haviour. If a parametrised system model can be contrived that is polynomial in form 
with the parameters appearing as coefficients, and the data can be transformed to repre- 
sent points on this polynomial, then any polynomial approximation method (including 
ours) may be used to estimate the parameters. 

In a case where several approximations exist, there is the question of how the char- 
acteristics of the approximation relate to the adequacy of the resulting system model. 
More specifically, it may be necessary to minimise some domain specific cost function 
that includes more than the error in the approximation and the contribution of higher 
order polynomial terms. By limiting the scope of this paper to the details of the CVB 
approximation algorithm per se, we leave such domain specific questions open for fu- 
ture discussion. 
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